 // ----------------------------------------------------CPS_taxsim_indexcalc.do
args reps indata cleandata moop_sims wealth_sims outdata

clear
local dsep=c(dirsep)
local reps=1

capture confirm file $ddCPS`dsep'CPS_taxsim_dataforfinalindex.dta BREAK
if _rc!=0 {

	use $ddCPS`dsep'CPS_forimputeandindex.dta, clear

	// Unless otherwise specified, year will be ref year
	rename yoc year
	
	// Bring in the MOOP data
	sort hid year
	cap drop _merge
	merge m:1 hid year using $ddCPS/CPS_mimp
	tab _merge
	tab year _merge
	drop _merge

	// Limit observations to matched sample
	keep if match==3

	/*
	cap assert m_sta1!=.
	if _rc!=0 {
		di "OBSERVATIONS NOT SUFFICIENTLY LIMITED IN THE IMPUTATION FILES"
		keep if m_sta1!=.
	}
	*/

	// Bring in the tax data!
	cap drop _merge 
	merge 1:1 yof h_seq ppos using CPS_fe_taxes
	assert _merge>1 
	tab yof _merge // there are lots of folks for whom we did taxsim who are not in the matched sample
	keep if _merge==3 
	
	// Keep only obs with non-zero weights
	drop if adj_pwgt==0|adj_pwgt==.
	
	// Bring in the wealth data
	sort hid year
	cap drop _merge
	merge m:1 hid using $ddCPS/CPS_wimp
	tab year _merge

	confirm var liq_`reps'
	confirm var m_sta`reps'
	confirm var m_dyn`reps'

	// Calculate income quintiles
		cap drop Y
		cap drop Ye
		qui gen Y=hinctot_tc if yof<1988
		qui replace Y=htotval_tc if yof>1987
		assert Y!=.
		qui gen Ye=Y/equiv
		// Pre drop
		xtset pid year
		qui gen l1_Ye=l1.Ye
		qui egen inc5=xtile(l1_Ye), by(year) weights(adj_pwgt) nq(5)
		label var inc5 "Pre-drop income quintiles"
		// Post drop
		qui egen Yeinc5=xtile(Ye), by(year) weights(adj_pwgt) nq(5)
		label var Yeinc5 "Y2 income quintiles"
		// Average
		bys pid: egen permYe=mean(Ye)
		qui egen permYeinc5=xtile(permYe), by(year) weights(adj_pwgt) nq(5)

	save $ddCPS`dsep'CPS_taxsim_dataforfinalindex.dta, replace
}



forvalues r=1(1)`reps' {
	quietly {
		timer clear
		timer on `r'
		use $ddCPS/CPS_taxsim_dataforfinalindex.dta, clear
		
		cap drop FW
		gen FW=liq_`r'

		cap drop tottax
		gen tottax = .5*hfica_bill_r`r' + hfiitax_r`r' + hsiitax_r`r' 
		// this fica var should be better than the calculated one from taxsim

		// Adjust for inflation
		cap drop cpi
		getcpi, year(year) gen(cpi)
		qui summ cpi if year==2010, mean
		replace cpi=cpi/`r(mean)'
		foreach var of varlist Y* FW tottax {
			replace `var'=`var'/cpi
		}
		
		// Topcode variables (temporary - done more sophisticated upstream)
		/*
		cap drop ntc_Y
		cap drop ntc_Ye
		rename Y ntc_Y
		rename Ye ntc_Ye
		qui gen Y=.
		qui gen Ye=.
		qui levelsof year, local(yrs)
		foreach y of local yrs {
			foreach var of varlist Y Ye {
				_pctile ntc_`var' if year==`y' [aw=adj_pwgt], p(98)
				qui replace `var'=min(ntc_`var',r(r1)) if year==`y'
			}
		} 
		*/


		// Calculate debt service
		local rn=.03
		cap drop DFW
		qui gen DFW=0 if FW>=0&FW<.
		qui replace DFW=FW*[`rn'/(1-((1+`rn')^-5))] if FW<0
		replace DFW=. if FW==.
		label var DFW "Debt service, r=`rn', time=5"

		// Code up income variables
		cap drop Yd
		qui gen Yd=(Y+DFW)/equiv
		cap drop tY 
		qui gen tY = Y-tottax
		cap drop tYe 
		qui gen tYe = (Y-tottax)/equiv
		cap drop tYd
		qui gen tYd=(Y+DFW-tottax)/equiv

		// Ye Adjusted for MOOP
		cap drop Ymei
		cap drop Ymeid
		qui gen Ymei=(Y*(1-m_sta`r'))/equiv
		qui gen Ymeid=(Y*(1-m_dyn`r'))/equiv
		cap drop Ymeis
		cap drop Ymeids
		rename Ymei Ymeis
		rename Ymeid Ymeids
		// MOOP and taxes
		cap drop tYmeis
		cap drop tYmeids
		qui gen tYmeis=(Y*(1-m_sta`r')-tottax)/equiv
		qui gen tYmeids=((Y*1-m_dyn`r')-tottax)/equiv
		// Ye adjusted for MOOP and debt service
		cap drop ymd0
		cap drop ymd1
		qui gen ymd0=Ymeis+(DFW/equiv) if Ymeis<.&DFW<.
		qui gen ymd1=Ymeids+(DFW/equiv) if Ymeids<.&DFW<.
		// Ye adjusted for MOOP, debt service, AND TAXES
		cap drop tymd0
		qui gen tymd0=Ymeis+((DFW-tottax)/equiv) if Ymeis<.&DFW<.&tottax<.
		cap drop tymd1
		qui gen tymd1=Ymeids+((DFW-tottax)/equiv) if Ymeids<.&DFW<.&tottax<.

		// Bottom code all income variables at 1
		foreach var of varlist Y* ym* tY* tym* {
			cap drop org_`var'
			gen org_`var'=`var'
			replace `var'=1 if `var'<1
		}

		// Calculate first differences
		tsset pid year
		foreach var of varlist Y Ye Yd tY tYe tYd {
			cap drop d`var'
			qui gen d`var'=d1.`var'/abs(l1.`var') if `var'<.&l1.`var'<.
		}
			// Pct change in income (w/ MOOP)
			cap drop dYmeicds
			qui gen dYmeicds=(Ymeids-l1.Ymeis)/abs(l1.Ymeis) if Ymeids<. & l1.Ymeis<.
			cap drop dtYmeicds
			qui gen dtYmeicds=(tYmeids-l1.tYmeis)/abs(l1.tYmeis) if Ymeids<. & l1.Ymeis<.
			// MOOP and debt service
			cap drop dymd
			qui gen dymd=(ymd1-l1.ymd0)/abs(l1.ymd0) if ymd1<.&l1.ymd0<.
			cap drop dtymd
			qui gen dtymd=(tymd1-l1.tymd0)/abs(l1.tymd0) if tymd1<.&l1.tymd0<.
			cap drop dYmd0
			cap drop dYmd1
			cap drop dYmd
			renpfix dy dY
			renpfix dty dtY 
			
		// Calculate ESIs
		foreach d in 4 2 10{
			local f`d'=-1/`d'
			foreach n of varlist dY* dtY* {
				loc v ""
				if substr("`n'",1,2)=="dY" loc v=subinstr("`n'","dY","Y",1)
				if substr("`n'",1,3)=="dtY" loc v=subinstr("`n'","dtY","tY",1)
				cap drop `v'_1o`d'_r`r'
				qui gen `v'_1o`d'_r`r' =(d`v'<=`f`d'') if d`v'<.
				label var `v'_1o`d'_r`r' "ESI 1/`d' for measure `v', rep `r'"
			}
		}
/*
		tabstat Ye_1o4 Ymeicds_1o4 Ymd_1o4 ///
			if !mi(Ye_1o4, Ymeicds_1o4, Ymd_1o4) [aw=adj_pwgt], by(year) f(%4.3f) s(mean)
*/			
		tabstat tYe_1o4 tYmeicds_1o4 tYmd_1o4 Ymd_1o4 ///
			if !mi(tYe_1o4, tYmeicds_1o4, tYmd_1o4, Ymd_1o4) [aw=adj_pwgt], by(year) f(%4.3f) s(mean)

		// Calculate drop tertiles and decide who can buffer
		foreach v of varlist Ymd_1o?_r`r' {
			local v1=subinstr("`v'","_r`r'","",.)
			cap drop tert_`v1'
			xtile tert_`v1'=dYmd if `v'==1 [aw=adj_pwgt], nq(3)
			replace tert_`v1'=4-tert_`v1'
			* tabstat dYmd [aw=adj_pwgt], by(tert_`v') s(min mean p50 max) f(%4.3f)
		}
		cap drop tert_1o?
		renpfix tert_Ymd_ tert_

		// Get head agegroup
		cap rename age_y age
		cap drop hagegroup
		gen hagegroup=.
		replace hagegroup=1 if hage<35
		replace hagegroup=2 if hage>=35 & hage<=44
		replace hagegroup=3 if hage>=45 & hage<=64
		replace hagegroup=4 if hage>=65 & hage<.
		label define hagegroup 1 "<35" 2 "35-44" 3 "45-64" 4 "65+", modify
		label val hagegroup hagegroup
/*
		// Income quintiles - MOVED
			// Pre drop
			xtset pid year
			qui gen l1_Ye=l1.Ye
			qui egen inc5=xtile(l1_Ye), by(year) weights(adj_pwgt) nq(5)
			label var inc5 "Pre-drop income quintiles"
			// Post drop
			qui egen Yeinc5=xtile(Ye), by(year) weights(adj_pwgt) nq(5)
			label var Yeinc5 "Y2 income quintiles"
			// Average
			bys pid: egen permYe=mean(Ye)
			qui egen permYeinc5=xtile(permYe), by(year) weights(adj_pwgt) nq(5)
*/
		// Bring in the buffers
		rename agegroup org_agegroup
		rename hagegroup agegroup
		sort tert_1o4 inc5 agegroup
		cap drop _merge
		merge tert_1o4 inc5 agegroup using ${ddCPS}/buffers_tert_1o4.dta, uniqusing
		* merge tert_1o4 inc5 agegroup using /home/shared/rehm/psid68-09/recovery/buffers_tert_1o4.dta // Old buffers
		rename agegroup hagegroup 

		rename org_agegroup agegroup 
		tsset pid year
		cap drop rwealth
		qui gen rwealth=(FW/l1.Ye)/equiv if l1.Ye<.
		qui summ rwealth [aw=adj_pwgt], d
		replace rwealth=r(p99) if rwealth>r(p99) & rwealth<.
		replace rwealth=r(p1) if rwealth<r(p1)
		capture rename buffer4_1o4 buffer6_1o4 
		local v "buffer6_1o4"
		cap drop wealthy_`v'
		gen wealthy_`v'=(rwealth>=`v') if rwealth<. & `v'<.
		label var wealthy_`v' "Dummy for being able to buffer inc drop"
		// Apply the wealth exclusion
		foreach v of varlist *Y*_1o4_r`r' {
			local j 6
			local n=substr("`v'",-1,1) // we don't do this anymore - only 1/4 counts
			cap drop w`j'`v'
			qui gen w`j'`v'=`v' if rwealth<.
			qui replace w`j'`v'=0 if w`j'`v'==1&wealthy_buffer`j'_1o4==1
		}

		// Next, exclude retirement entry
		foreach v of varlist *_1o?_r`r' {
			cap drop r`v'
			qui gen r`v'=`v'
			* qui replace r`v'=0 if r`v'==1&(heret==1|spret==1)
			qui replace r`v'=0 if r`v'==1&hesp_ret==1
			* qui replace r`v'=. if hesp_ret==1
		}



		// Create some vars to collapse with
		cap drop period
		qui gen period=.
		qui replace period=0 if year<=1996
		qui replace period=1 if year>1996&year<=2007
		qui replace period=2 if year>2007
		label define period 0 "1986-1996" 1 "1997-2007" 2 "2008-2010", replace
		label val period period

		
		confirm var adults
		confirm var kids
		cap drop HH_type
		gen HH_type=.
		replace HH_type=1 if (adults==1 & kids==0)
		replace HH_type=2 if (adults==1 & kids>0)
		replace HH_type=3 if (adults>1 & kids==0)
		replace HH_type=4 if (adults>1 & kids>0)
		label define HH_type 1 "single individual" 2 "single parent" 3 "several adults w/o kids" 4 "several adults w/kids", modify
		label val HH_type HH_type
		
		
		cap drop sample_r`r'
		qui gen sample_r`r'=1 if !mi(rYe_1o4_r`r', rYmeicds_1o4_r`r', rYmd_1o4_r`r', rw6Ymd_1o4_r`r')
/*
		tab sample, mi
		tabstat 			/// 
			rw6Ye_1o4 		/// income+buffer(exclusion)
			rw6Ymeicds_1o4 	/// income+moop+buffer
			rw6Ymd_1o4 		/// income+moop+debt+buffer
			if sample==1 [aw=adj_pwgt], by(year) f(%4.3f) s(mean) 

		// Median drop
		tabstat dYmd if rw6Ymd_1o4==1 [aw=adj_pwgt], s(median) f(%5.4f) by(year)
		tabstat dYmd if rw6Ymd_1o4==1 [aw=adj_pwgt], s(median) f(%5.4f) by(period)
*/
STOP
	keep yof h_seq ppos *_r`r' 
	sort yof h_seq ppos
	tempfile t`r'
	save `t`r'', replace
	
	timer off `r'
	timer list `r'
	}
	di "`r' in `r(t`r')' seconds"
}
forvalues r=1/`reps' {
	if `r'==1 use `t1', clear
	else {
		cap drop _merge
		merge 1:1 yof h_seq ppos using `t`r''
	}
}
		save $ddCPS/CPS_taxsim_esidata.dta, replace

		//------------------------------------------------------------
		// Collapse drops
		//------------------------------------------------------------
		cap rename heduc_y heduc
		keep *1o? dYmd adj_pwgt year period sample hid pid ///
			Ye permYe age agegroup heduc hrace sex permYeinc5 HH_type ///
			state
		keep if sample==1
		drop if rw6Ymd_1o4==.
/*		
		table agegroup period if sample==1 [aw=adj_pwgt],	///
			cont(mean rw6Ymd_1o4) f(%4.3f)
*/
		gen income = ""
		forvalues j=1/5{
			replace income="average HH Ye for xt5==`j'" if permYeinc5==`j'
		}
		// age
		drop age
		label define agegroup 0 "<17" 1 "18-34" 2 "35-44" 3 "45-64" 4 "65+", modify

		label val agegroup agegroup
		decode agegroup, gen(age)
		*/
		// edu
		* qui replace heduc=2 if heduc==3
		* qui replace heduc=4 if heduc==5
		label define heduc 1 "LT_HS" ///
						   2 "HS"	 ///
						   3 "S_COL" ///
						   4 "COLL"	 ///
						   5 "P_COL" 
		label val heduc heduc
		cap drop educ
		decode heduc, gen(educ)
		// Race
		label define hrace  1 "White NH" ///
							2 "Black NH" ///
							3 "Other NH" ///
							4 "Hispanic"
		label val hrace hrace
		decode hrace, gen(race)
		// Gender
		qui gen gender="female" if sex==2
		qui replace gender="male" if sex==1
		// HH type
		decode HH_type, gen(famtype)
		foreach v of varlist *Y*_1o? {
			gen obs`v'=`v'
		}
		// State
		#d ;
		label values state  state;
		label define state
			11          "ME"                         
			12          "NH"                 
			13          "VT"                       
			14          "MA"                 
			15          "RI"                  
			16          "CT"                   
			21          "NY"                      
			22          "NJ"                    
			23          "PA"                  
			31          "OH"                          
			32          "IN"                       
			33          "IL"                      
			34          "MI"                      
			35          "WI"                     
			41          "MN"                     
			42          "IA"                          
			43          "MO"                      
			44          "ND"                  
			45          "SD"                  
			46          "NE"                      
			47          "KS"                        
			51          "DE"                      
			52          "MD"                      
			53          "DC"          
			54          "VA"                      
			55          "WV"                 
			56          "NC"                
			57          "SC"                
			58          "GA"                       
			59          "FL"                       
			61          "KY"                      
			62          "TN"                     
			63          "AL"                       
			64          "MS"                   
			71          "AR"                      
			72          "LA"                     
			73          "OK"                      
			74          "TX"                         
			81          "MT"                       
			82          "ID"                         
			83          "WY"                       
			84          "CO"                      
			85          "NM"                    
			86          "AZ"                       
			87          "UT"                          
			88          "NV"                        
			91          "WA"                    
			92          "OR"                        
			93          "CA"                    
			94          "AK"                        
			95          "HI"                        
		;
		#d cr
		cap drop temp
		cap drop region
		qui gen temp = floor(state/10)
		qui gen region=.
		qui replace region=1 if temp<3&temp>0
		qui replace region=2 if temp>2&temp<5
		qui replace region=3 if temp>4&temp<8
		qui replace region=4 if temp>7&temp<10
		assert region>0&region<5
		label define region 1 "North East" 2 "North Central" 3 "South" 4 "West"
		label val region region 
		rename temp division
		label define div 1 "New England" 2 "Mid-Atlantic" 3 "E North Central" 4 "W North Central" ///
			5 "South Atlantic" 6 "East South" 7 "West South" 8 "Mountain" 9 "Pacific"
		label val division div
		cap drop tempd
		cap drop tempr
		rename division tempd
		rename region tempr
		rename state hg_st60
		decode hg_st60, gen(state)
		decode tempd, gen(division)
		decode tempr, gen(region)

		order year adj_pwgt age educ famtype gender income race state region division *Y*_1o? obs*
		// COLLAPSE
		gen all="all"
		local c "all age educ race gender income famtype state region division"
		foreach v of local c {
			di "variable: `v'"
			preserve
				collapse (mean) r*Y*_* (count) obs* [aw=adj_pwgt], by(year `v')
				gen GROUP="`v'"
				save "$ddCPS`dsep'CPS_results_`v'.dta", replace
			restore
		}
		*/
		use "$ddCPS/CPS_results_all.dta", clear
		qui gen group=all
		foreach f in age educ race gender income famtype state region division {
			append using "$ddCPS/CPS_results_`f'.dta"
			qui replace GROUP="`f'" if GROUP==""
			replace group=`f' if group==""
			drop if group==""
			drop `f'
		}
		drop obs*
		drop all		
		order GROUP group year

		// Add in period
		cap drop period
		qui gen period=.
		qui replace period=0 if year<=1996
		qui replace period=1 if year>1996&year<=2007
		qui replace period=2 if year>2007
		label define period 0 "1986-1996" 1 "1997-2007" 2 "2008-2010"
		label val period period
		// State period
		qui gen stperiod=.
		qui replace stperiod=0 if year>2004&year<2008
		qui replace stperiod=1 if year>2007
		label define stperiod 0 "2005-2007" 1 "2008-2010"
		label val stperiod stperiod

		// Collapse period drops
		tostring year, gen(YEAR)

		preserve
			collapse (mean) r*Y*_1o4, by(GROUP group period)
			decode period, gen(YEAR)
			cap drop period
			qui gen period=1
			tempfile t1
			save `t1', replace
		restore

		preserve
			collapse (mean) r*Y*_1o4, by(GROUP group stperiod)
			decode stperiod, gen(YEAR)
			cap drop stperiod
			qui gen period=2
			tempfile t2
			save `t2', replace
		restore

		cap drop period
		qui gen period=0

		append using `t1'
		append using `t2'
		drop year
		rename YEAR year
		* qui drop if period==2&GROUP!="state"
		* qui drop if GROUP=="state"&period!=2
		qui drop if year==""
		compress
		sort GROUP group period year
		save $ddCPS/CPS_ESIresults_r`r'.dta, replace
		timer off `r'
	}
	di "`r'"
	timer list `r'
}
exit
outsheet using $oCPS/Index_final.csv, name comma replace
exit
